Eliminación de ruido en imágenes con wavelets.

Análisis de señales

<<<<<<< HEAD ======= >>>>>>> Alejandra

1 Introducción

En el campo del análisis de señales, uno de los problemas más relevantes es la eliminación de ruido en imágenes digitales. Las señales, que representan variaciones de magnitudes en el espacio y/o tiempo, a menudo están contaminadas por ruido debido a interferencias en los procesos de captura o transmisión. Este ruido, que puede existir de formas muy diversas como ruido gaussiano, artefactos asociados a frecuencias altas o bajas, o efectos multiplicativos, complica la extracción de características significativas.

El uso de wavelets es una herramienta clave para abordar este problema, el cual permite descomponer señales en sus componentes espaciales y frecuenciales de forma eficiente. Este método ofrece la posibilidad de adaptarse a las características particulares de cada tipo de ruido.

En este trabajo se propone explorar la capacidad de los wavelets para la eliminación de ruido en imágenes digitales. Para ello, se generarán y analizarán diferentes tipos de ruido artificial, evaluando su complejidad para su atenuación y determinando los parámetros más eficaces de los wavelets.

2 Fundamento teórico: eliminación de ruido con wavelets

Los wavelets son funciones matemáticas que permiten descomponer una señal en componentes de diferentes escalas, lo que resulta útil para identificar y procesar características específicas. Este enfoque es especialmente relevante en la eliminación de ruido, donde las frecuencias indeseadas pueden ser separadas y reducidas sin afectar significativamente las características principales de la señal original. A diferencia de la Transformada de Fourier, que opera globalmente y no ofrece información sobre la localización temporal de los eventos, los wavelets permiten un análisis localizado, facilitando la identificación de patrones y anomalías en los datos.

El proceso de reducción de ruido con wavelets generalmente incluye tres etapas principales: la descomposición de la señal utilizando una wavelet madre, la modificación de los coeficientes wavelet mediante técnicas de umbral, y la reconstrucción de la señal. La selección de la wavelet madre adecuada y los parámetros de umbral son aspectos críticos que deben adaptarse a las características específicas del ruido y la señal.

Además, los wavelets ofrecen un enfoque multi-resolución, permitiendo una representación detallada de los componentes de alta frecuencia, asociados frecuentemente con el ruido, mientras preservan las estructuras globales de baja frecuencia. Esta característica hace que los wavelets sean especialmente útiles en aplicaciones donde la precisión y la integridad de los datos son esenciales, como en imágenes médicas, procesamiento de audio o análisis de datos científicos.

En este trabajo, se emplearán wavelets como herramienta principal para eliminar diferentes tipos de ruido sintético en imágenes.

3 Funciones de programacion empleadas

4 Desarrollo y resultados

5 Fundamento teórico: eliminación de ruido con wavelets

6 Funciones de programacion empleadas

Eliminación de ruido con el algoritmo de Mallat

Empleamos la funicón imwd() del paquete wavethresh. Está función realiza una transformada discreta wavelet de acuerdo con el algoritmo de Mallat.

• imwd(image, filter.number=10, family="DaubLeAsymm", type="wavelet",...).

El argumento image de la función debe ser una matriz cuadrada cuya dimensión sea potencia de dos.filter.number elige la suavidad de la wavelet a emplear, siendo por defecto 10. family indica la familia de wavelets a emplear (“DaubExPhase” ó “DaubLeAsymm”).Para tratar las fronteras, mantendremos el parámetro bc = "periodic" por defecto.

Para la eliminación de ruido aplicaremos la función threshold() al objeto que devuelve la función imwd().

• threshold(imwd, levels = 3:(nlevelsWT(imwd) - 1), type = "hard", policy = "universal", by.level = FALSE, value = 0, return.threshold = FALSE, compression = TRUE, Q = 0.05, ...)

levels es el número de niveles a los cuáles deseamos aplicar un umbral, mientras que type indica si queremos un umbral más suave (“soft”) ó más fuerte (“hard”). El parámetro policy selecciona la técnica para elegir umbral. by.level = FALSE significa que se aplica un umbral global a todos los niveles indicados, mientras que by.level = TRUE calcular threshold para cada nivel por seaparado. El parámetro value, el valor del umbral, se usará si elegimos técnicas manuales en policy. Si queremos obtener el valor umbral aplicado, pondremos return.threshold = TRUE. Finalmente, dejaremos el parámetro compression = TRUE por defecto, para obterner un objeto más pequeño.

Durante el desarrollo del trabajo variaremos estos parámetros para observar su efecto en la eliminación de diversos tipos de ruido y encontar aquellos que generen un mejor resultado en cada caso. Una vez realizado el thresholding, emplearemos la función imwr para realizar la transformda wavelet inversa y poder visualizar la imagen tras la eliminación de ruido.

7 Desarrollo y resultados

Antes de comenzar, instalamos y/o cargamos todos los paquetes requeridos.

if (!require("pacman")) install.packages("pacman")
library(pacman)
p_load(imager, wavethresh, ggplot2, dplyr, SpatialPack, waveslim, EBImage, stringr, jpeg, abind, magick)

Comenzamos cargando y visualizando las fotografías a emplear.

images_path <- list.files("./fotos", full.names = TRUE)

nombres_images <- str_remove_all(string = str_remove_all(string = images_path, pattern = "./fotos/"), pattern = "\\.JPG|\\.jpg")

images <- lapply(images_path, readJPEG) # Cargamos las imágenes
names(images) <- nombres_images

# Rotamos algunas de las fotos para una visualización más uniforme
fotos_a_girar <- c("1", "2", "3", "4")
images_rotadas <- lapply(images[fotos_a_girar], aperm, perm = c(2, 1, 3))


for (i in fotos_a_girar) {
  images_rotadas[[i]] < images_rotadas[[i]][dim(images_rotadas[[i]])[1]:1, , ]
}

images[fotos_a_girar] <- images_rotadas
rm(images_rotadas)
rm(fotos_a_girar)

# Visualizamos las imagenes originales (comento para que no tarde en ejecutar)

#par(mfrow = c(2, 3), mar = c(1, 1, 1, 1))

# for (img in nombres_images) {
#   display(Image(images[[img]], colormode = "Color"), method = "r")
#   title(paste('Imagen', img))
# }

7.1 Inclusión de ruido sintético en las imágenes

<<<<<<< HEAD ======= >>>>>>> Alejandra

# Definición de tipos de ruido

NOISE_TYPES <- list(
  gaussian = list(
    generator = function(channel, params) {
      # Desviación estándar del ruido con un valor predeterminado
      noise_std_dev <- params$std_dev %||% 0.5

      # Generación de ruido gaussiano
      noise <- array(
        rnorm(length(channel), mean = 0, sd = noise_std_dev),
        dim = dim(channel)
      )

      # Asegurar que los valores estén entre 0 y 1
      pmax(0, pmin(1, channel + noise))
    }
  ),
  sinusoidal_high = list(
    generator = function(channel, params) {
      # Frecuencia y amplitud del ruido sinusoidal de alta frecuencia
      frequency <- params$frequency %||% 25
      amplitude <- params$amplitude %||% 0.2

      # Generación de ruido sinusoidal
      height <- dim(channel)[1]
      width <- dim(channel)[2]
      x <- seq(0, 2 * pi, length.out = width)
      y <- seq(0, 2 * pi, length.out = height)
      noise_grid <- outer(sin(x * frequency), sin(y * frequency))

      # Aplicar el ruido
      noise <- array(noise_grid * amplitude, dim = dim(channel))
      pmax(0, pmin(1, channel + noise))
    }
  ),
  sinusoidal_low = list(
    generator = function(channel, params) {
      # Frecuencia y amplitud del ruido sinusoidal de baja frecuencia
      frequency <- params$frequency %||% 2
      amplitude <- params$amplitude %||% 0.2

      # Generación de ruido sinusoidal
      height <- dim(channel)[1]
      width <- dim(channel)[2]
      x <- seq(0, 2 * pi, length.out = width)
      y <- seq(0, 2 * pi, length.out = height)
      noise_grid <- outer(sin(x * frequency), sin(y * frequency))

      # Aplicar el ruido
      noise <- array(noise_grid * amplitude, dim = dim(channel))
      pmax(0, pmin(1, channel + noise))
    }
  ),
  salt_pepper = list(
    generator = function(channel, params) {
      # Proporción de píxeles afectados por el ruido de sal y pimienta
      epsilon <- params$epsilon %||% 0.2

      # Generación de ruido
      noise <- matrix(sample(c(0, 1, NA), length(channel), replace = TRUE, prob = c(epsilon / 2, epsilon / 2, 1 - epsilon)),
        nrow = dim(channel)[1], ncol = dim(channel)[2]
      )
      channel[!is.na(noise)] <- noise[!is.na(noise)]
      channel
    }
  ),
  gamma = list(
    generator = function(channel, params) {
      # Ruido multiplicativo gamma con parámetro de dispersión
      looks <- params$looks %||% 2
      noise <- array(rgamma(length(channel), shape = looks, scale = 1 / looks), dim = dim(channel))
      pmax(0, pmin(1, channel * noise))
    }
  ),
  uniform_multiplicative = list(
    generator = function(channel, params) {
      # Ruido multiplicativo uniforme
      looks <- params$looks %||% 2
      noise_channel <- SpatialPack::imnoise(
        img = channel,
        type = "speckle",
        looks = looks
      )
      pmax(0, pmin(1, noise_channel))
    }
  )
<<<<<<< HEAD
)

# Función para añadir ruido a una imagen
add_noise_to_image <- function(image_name, noise_type, noise_params = list(), plot = FALSE) {
=======
)

# Función para añadir ruido a una imagen
add_noise_to_image <- function(image_name, noise_type, noise_params = list(),plot=FALSE) {
>>>>>>> Alejandra
  # Verificar si la imagen existe en la lista
  if (!image_name %in% names(images)) {
    stop("La imagen con este nombre no se encuentra en la lista 'images'")
  }

  # Verificar el tipo de ruido
  if (!noise_type %in% names(NOISE_TYPES)) {
    stop(
      "El tipo de ruido es desconocido. Tipos disponibles: ",
      paste(names(NOISE_TYPES), collapse = ", ")
    )
  }

  # Obtener la imagen original de la lista
  original_image <- images[[image_name]]

  # Convertir la imagen a un array si es necesario
  image_array <- as.array(original_image)

  # Aplicar ruido a cada canal
  noisy_channels <- lapply(1:3, function(i) {
    channel <- image_array[, , i]
    NOISE_TYPES[[noise_type]]$generator(channel, noise_params)
  })

  # Crear la imagen con ruido
  noisy_image_array <- array(
    unlist(noisy_channels),
    dim = dim(image_array)
  )

<<<<<<< HEAD
  # Visualizar si se ha indicado
=======
 # Visualizar si se ha indicado
>>>>>>> Alejandra
  if (plot == TRUE){
  layout(matrix(1:2, 1, 2))
  plot(Image(original_image, colormode = "Color"))
  title("Original")
  plot(Image((noisy_image_array), colormode = "Color"))
  title(paste("Ruido:", noise_type))}
  
  return(noisy_image_array)
}
<<<<<<< HEAD


# Aplicar los diferentes tipos de ruido a cada imagen
add_noise_to_image("1", "gaussian", list(std_dev = 0.3))

 [ reached getOption("max.print") -- omitted 3 matrix slice(s) ]

add_noise_to_image("2", "sinusoidal_high", list(frequency = 25, amplitude = 0.2))

 [ reached getOption("max.print") -- omitted 3 matrix slice(s) ]

add_noise_to_image("3", "sinusoidal_low", list(frequency = 2, amplitude = 0.2))

 [ reached getOption("max.print") -- omitted 3 matrix slice(s) ]

add_noise_to_image("4", "salt_pepper", list(epsilon = 0.1))

 [ reached getOption("max.print") -- omitted 3 matrix slice(s) ]

add_noise_to_image("5", "gamma", list(looks = 2))

 [ reached getOption("max.print") -- omitted 3 matrix slice(s) ]

add_noise_to_image("5", "uniform_multiplicative", list(looks = 2))

 [ reached getOption("max.print") -- omitted 3 matrix slice(s) ]

# Representación de las dos familias de wavelets empleadas por la función imwd.

# Generación de un filtro con los coeficientes correspondientes de cada wavelet
wave_coeffs <- filter.select(filter.number = 10, family = "DaubLeAsymm")
wave_coeffs_2 <- filter.select(filter.number = 10, family = "DaubExPhase")

# Eje temporal
x <- seq_along(wave_coeffs$H)

# # Visualización de las wavelets
par(mfrow = c(1, 2))

plot(x, wave_coeffs$H, type = "l", lwd = 2,
     main = "Wavelet DaubLeAsymm",
     xlab = "x", ylab = "Amplitud")
points(x, wave_coeffs$H, pch = 19)

plot(x, wave_coeffs_2$H, type = "l", lwd = 2,
     main = "Wavelet DaubExPhase",
     xlab = "x", ylab = "Amplitud")
points(x, wave_coeffs_2$H, pch = 19)

rm(wave_coeffs)
rm(wave_coeffs_2)
rm(x)

=======

# Aplicar los diferentes tipos de ruido a cada imagen
#add_noise_to_image("1", "gaussian", list(std_dev = 0.3))
#add_noise_to_image("2", "sinusoidal_high", list(frequency = 25, amplitude = 0.2))
#add_noise_to_image("3", "sinusoidal_low", list(frequency = 2, amplitude = 0.2))
#add_noise_to_image("4", "salt_pepper", list(epsilon = 0.1))
#add_noise_to_image("5", "gamma", list(looks = 2))
#add_noise_to_image("5", "uniform_multiplicative", list(looks = 2))

>>>>>>> Alejandra

<<<<<<< HEAD

7.2 Función imwd

El primer método que emplearemos para la eliminación de ruido es el algoritmo de Mallat a través de la función imwd. Antes de poder aplicar esta función, necessitames realizar un pre-procesamiento a las imágenes, ya que este algoritmo necesita una matriz cuadrada cuyas dimensiones sean potencia de dos. Se han escogido 2 maneras distintas para obtener imágenes con el tamaño adecuado. Por un lado, redimensionaremos las imágenes con la función resize, lo que podría conllevar problemas de distorsión si las imágenes estaban lejos de tener dimensiones cuadradas. Por ello, también vamos a emplear otra técnica y rellenaremos las matrices de las imágnes con valores de 0 hasta alcanzar las dimensiones adecuadas. Observamos los problemas que surguen en cada caso y trataremos de solucionarlos de distintas maneras.

7.2.1 Redimensionando las imágenes

Como ya se ha visto para poder aplicar la función imwd()es necesario partir de una matriz cuadrada cuyas dimensiones sean potencia de dos. Por ello, en primer lugar creamos una función resize_imwd() tal que dada una foto busca la submatriz cuadrada y potencia de dos más grande posible y a continuacón redimensiona la imagen a dicha submatriz cuadrada.

# Pre-procesamiento al aplicado de función imwd (algoritmo de Mallat).
# Redimensionamiento de la imagen 

resize_imwd<- function(foto){
  
  img <- as.cimg(foto) # Pasar a formato Imager para aplicar función resize
   
  # Dimensiones de la foto
  dim_foto <- dim(foto)
  filas <- dim_foto[1]
  columnas <- dim_foto[2]
  
  lado_minimo <- min(filas, columnas)  # Tamaño submatriz cuadrada mas grande
  lado_potencia2 <- 2^floor(log2(lado_minimo)) # Tamaño submatriz cuadrada potencia de dos mas grande
  
  foto_resized <-resize(foto, w = lado_potencia2, h = lado_potencia2) # Redimensionado de la imagen

return(foto_resized)
}

Por otro lado, creamos una función para el post-procesamiento de las imágenes tras la eliminición de ruido. Queremos devolverlas a su tamaño original con el objetivo de comparar con las imágenes iniciales.

# Post-procesamiento de la imagen:
# Redimensionamiento de la imagen a su tamaño original.

resize_imwd_to_original<- function(image_redimensionada, nombre_foto){
  
  #img <- as.cimg(image_redimensionada)
  foto <- images[[nombre_foto]]
  
  # Dimensiones de la foto
  dim_foto <- dim(foto)
  filas <- dim_foto[2]
  columnas <- dim_foto[1]
  
  # Redimensionado de la imagen
  foto_resized <-EBImage::resize(image_redimensionada, w = columnas, h = filas) 
 

return(foto_resized)
}

Comenzamos generando una función procesar_imagen_wavelet con parámetros foto, tipo y policy. Esta función realiza en primer lugar la transformada wavelet a cada uno de los tres canales de una imagen. A continuación, se realiza el thresholding con la función threshold, pudiendo variar de el tipo de “hard” a “soft” y el parámetro policy (modificando adecuadamente los parámetros necesarios en la función threshold en este último caso). Una vez realizada la eliminación de ruido, se aplica la trasnformada wavelet inversa imwr para por último reconstruir la imagen a partir de los tres canales.

procesar_imagen_wavelet <- function(foto, tipo = "hard", policy = "universal") {
  # 1. Realizamos la transformada wavelet a cada canal
  lwd <- lapply(1:3, function(canal) {
    imwd(foto[,,canal])  
  })
  
  # 2. Aplicamos el umbral a los coeficientes de la transformada wavelet
  lwd_threshold <- lapply(lwd, function(canal_wd) {
    niveles <- canal_wd$nlevels
    wavethresh::threshold(canal_wd, levels = 3:(niveles-1), type = tipo, policy = policy,by_level=TRUE,compression=FALSE)
  })
  # 3. Aplicamos la transformada wavelet inversa a cada canal umbralizado
  ilwd <- lapply(lwd_threshold, function(canal_umbralizado) {
    wavethresh::imwr(canal_umbralizado)  # Transformada wavelet inversa
  })
  
  # 4. Reconstruir la imagen combinando los tres canales procesados
  imagen_reconstruida <- abind::abind(ilwd[[1]], ilwd[[2]], ilwd[[3]], along = 3)
    imagen <- Image(imagen_reconstruida, colormode = 'Color')
  
  return(imagen)
}

Para comenzar el análisis, vamos a emplear dos imágenes muy parecidas (imágenes 4 y 5). Una de ellas tiene muy alta resolución mientras que la segunda cuenta con una calidad mucho menor. El objetivo es determinar si la resolución de la imagen afecta a la hora de eliminar ruido de esta. Vamos a probar con el primer tipo de ruido, ruido gaussiano.

# Añadimos ruido gaussiano a las imágenes con sd = 0.3
image_4_gaussian_noise <- add_noise_to_image("4", "gaussian", list(std_dev = 0.3))
image_5_gaussian_noise <- add_noise_to_image("5", "gaussian", list(std_dev = 0.3))

# Hacemos una lista con las imágenes con ruido a redimensionar
images_for_imwd_cut_1 <- list(image_4_gaussian_noise, image_5_gaussian_noise)
names(images_for_imwd_cut_1) <- c('4 Noise: gaussian', '5 Noise: gaussian')

# Eliminamos variables innecesarias
rm(image_4_gaussian_noise)
rm(image_5_gaussian_noise)

Aplicamos la función resize_imwd creada anteriormente para obtener una matriz con las dimensiones necesarias para aplicar la transformada wavelet.

images_recortadas <- lapply(images_for_imwd_cut_1, resize_imwd)

Warning: Assuming third dimension corresponds to colourWarning: Assuming third dimension corresponds to colour

Una vez tenemos las imágenes con ruido generadas y redimensionadas adecuadamente, podemos aplicar la función imwd a cada uno de los tres canales (R, G y B). Empleamos la función procesar_imagen_wavelet que devuelve las imágenes reconstruidas después del thresholding. Para ruido gaussiano y para comenzar, vamos a elegir los parámetros por defecto de la función para el thresholding.

images_sin_ruido_gaussiano <- lapply(images_recortadas, procesar_imagen_wavelet)

Finalmente, visualizamos los resultados. Primeramente, observamos las imágenes redimensionadas con ruido y la imagen obtenida tras el uso del método de thresholding para la eliminación de este. Observamos una principal diferencia entre ambas: la foto que contaba con menor resolución presenta también el peor resultado. Aunque el ruido haya sido eliminadom, sus bordes están más difuminados y tiene muy baja calidad.

par(mfrow = c(2, 2), cex = 0.5)

display(Image(images_recortadas[[1]], colormode = 'Color'), method='r')
title('Imagen 4 con ruido')

display(Image(images_sin_ruido_gaussiano[[1]], colormode = 'Color'), method='r')
title('Imagen 4 sin ruido')

display(Image(images_recortadas[[2]], colormode = 'Color'), method='r')
title('Imagen 5 con ruido')

display(Image(images_sin_ruido_gaussiano[[2]], colormode = 'Color'), method='r')
title('Imagen 5 sin ruido')

En segundo lugar, vamos a visualizar las imágenes originales y las imágenes sin ruido redimensionadas a su tamaño original, usando la función resize_imwd_to_original.

images_sin_ruido <- Map(resize_imwd_to_original, images_sin_ruido_gaussiano, c("4", "5") )

par(mfrow = c(2, 2), cex = 0.5)

display(Image(images[[4]], colormode = 'Color'), method='r')
title('Imagen 4 original')

display(Image(images_sin_ruido[[1]], colormode = 'Color'), method='r')
title('Imagen 4 sin ruido')

display(Image(images[[5]], colormode = 'Color'), method='r')
title('Imagen 5 original')

display(Image(images_sin_ruido[[2]], colormode = 'Color'), method='r')
title('Imagen 5 sin ruido')

Vemos que efectivamente, la imagen que originalemente contaba con un número de píxeles mucho menor, la imagen 5, presenta una gran distorsión tras la eliminación de ruido.

rm(images_for_imwd_cut_1,images_sin_ruido_gaussiano)

A continuación vamos a comprobar que es lo que ocurre cuando añadimos ruido sintético sinusoidal y si la frecuencia de este afecta al resultado de la eliminación de ruido. Trabajaremos con una única fotografía: la imagen 1.

# Añadimos ruido sinusoidal a las imágenes con sd = 0.3
image_1_sinosuidal_high <- add_noise_to_image("1", "sinusoidal_high", list(frequency = 50, amplitude = 0.3))
image_1_sinosuidal_low <- add_noise_to_image("2", "sinusoidal_low", list(frequency = 5, amplitude = 0.3))

# Hacemos una lista con las imágenes con ruido a redimensionar
images_for_imwd_cut_2 <- list(image_1_sinosuidal_high , image_1_sinosuidal_low )
names(images_for_imwd_cut_2) <- c('1 Noise: sinusoidal high', '1 Noise: sinusoidal low')

# Eliminamos variables innecesarias
rm(image_1_sinosuidal_high)
rm(image_1_sinosuidal_low)

images_recortadas <- lapply(images_for_imwd_cut_2, resize_imwd)

Warning: Assuming third dimension corresponds to colourWarning: Assuming third dimension corresponds to colour

images_sin_ruido_sinusoidal<- lapply(images_recortadas, procesar_imagen_wavelet)

Esta claro que los parámetros por defecto de la función threshold no son capaces de eliminar el ruido de tipo sinusoidal de la manera en que si lo era con el ruido de tipo Gaussiano, un tipo de ruido aleatorio, al contrario que el sinusoidal, que es una señal periódica.

par(mfrow = c(2, 2), cex = 0.5)

display(Image(images_recortadas[[1]], colormode = 'Color'), method='r')
title('Imagen 1 con ruido de frecuencia alta')

display(Image(images_sin_ruido_sinusoidal[[1]], colormode = 'Color'), method='r')
title('Imagen 1 sin ruido')

display(Image(images_recortadas[[2]], colormode = 'Color'), method='r')
title('Imagen 1 con ruido de frecuencia baja')

display(Image(images_sin_ruido_sinusoidal[[2]], colormode = 'Color'), method='r')
title('Imagen 1 sin ruido')

Vamos a variar parámetros de la función threshold para intentar quitar este ruido de manera más manual. En primer lugar, cambiamos el número de niveles al que aplicamos el umbral, para incluirlos a todos. Cambiamos policy a “manual” y variamos el valor del umbral de forma manual hasta encontrar uno que sea satisfactorio.

procesar_imagen_wavelet_sinusoidal <- function(foto, tipo = "hard", policy = "universal") {
  # 1. Realizamos la transformada wavelet a cada canal
  lwd <- lapply(1:3, function(canal) {
    imwd(foto[,,canal])  
  })
  
  # 2. Aplicamos el umbral a los coeficientes de la transformada wavelet
  lwd_threshold <- lapply(lwd, function(canal_wd) {
    niveles <- canal_wd$nlevels
    wavethresh::threshold(canal_wd, levels = 1:(niveles-1), type = tipo, policy = policy,by_level=TRUE,compression=FALSE, value = 4)
  })
  # 3. Aplicamos la transformada wavelet inversa a cada canal umbralizado
  ilwd <- lapply(lwd_threshold, function(canal_umbralizado) {
    wavethresh::imwr(canal_umbralizado)  # Transformada wavelet inversa
  })
  
  # 4. Reconstruir la imagen combinando los tres canales procesados
  imagen_reconstruida <- abind::abind(ilwd[[1]], ilwd[[2]], ilwd[[3]], along = 3)
    imagen <- Image(imagen_reconstruida, colormode = 'Color')
  
  return(imagen)
}

images_sin_ruido_sinusoidal<- lapply(images_recortadas, procesar_imagen_wavelet_sinusoidal, tipo ="soft", policy = "manual")

Con un valor de umbral de 4 y variando el tipo a “soft”, observamos que hemos conseguido eliminar el ruido sinusoidal de alta frecuencia, pero pagando un precio muy alto: los bordes de la imágen se disorsionan completamente y tenemos una muy baja resolución. Por otro lado, el ruido de frecuencia baja, aunque ha disminuido, claramente sigue presente en la imagen. El umbral necesario para eliminarlo con este método es tan alto que distorsionaría la imagen casi por completo.

par(mfrow = c(2, 2), cex = 0.5)

display(Image(images_recortadas[[1]], colormode = 'Color'), method='r')
title('Imagen 1 con ruido de frecuencia alta')

display(Image(images_sin_ruido_sinusoidal[[1]], colormode = 'Color'), method='r')
title('Imagen 1 sin ruido')

display(Image(images_recortadas[[2]], colormode = 'Color'), method='r')
title('Imagen 1 con ruido de frecuencia baja')

display(Image(images_sin_ruido_sinusoidal[[2]], colormode = 'Color'), method='r')
title('Imagen 1 sin ruido')

Probamos otros tipos de ruido: gamma, salt and pepper y ruido uniforme en las imágenes 2 y 3.

# Añadimos ruido a imágenes
image_2_salt_pepper <- add_noise_to_image("2", "salt_pepper", list(epsilon = 0.1))
image_3_gamma <- add_noise_to_image("3", "gamma", list(looks = 2))
image_3_uniform <- add_noise_to_image("3", "uniform_multiplicative", list(looks = 2))


# Hacemos una lista con las imágenes con ruido a recortar
images_for_imwd_cut_3 <- list( image_2_salt_pepper, image_3_gamma, image_3_uniform)
names(images_for_imwd_cut_3) <- c('2 Noise: salt and pepper', '3 Noise: gamma','3 Noise: uniform')

# Eliminamos variables innecesarias
rm(image_2_salt_pepper, image_3_gamma, image_3_uniform)
=======

4.2 Función imwd








Función para hacer cuadrada la imagen.
Para aplicar el algoritmo de Deformación Iterativa de Mallat (IMWD),
es necesario que la imagen tenga una forma cuadrada. Dado que muchas
imágenes no son cuadradas, es necesario convertirlas antes de aplicar el
algoritmo.
A continuación, se presenta una función de pre-procesamiento, que
ajusta cualquier imagen rectangular a un tamaño cuadrado, manteniendo
sus proporciones originales al agregar relleno (fondo blanco) si es
necesario. Esta transformación asegura que la imagen sea compatible con
el algoritmo IMWD.




hacer_cuadrada_potencia_2 <- function(imagen) {
  n_filas <- dim(imagen)[1]
  n_columnas <- dim(imagen)[2]
  
  nuevo_tamano <- max(n_filas, n_columnas)
  
  siguiente_potencia_2 <- 2^ceiling(log2(nuevo_tamano))
  
  imagen_cuadrada <- array(0, dim = c(siguiente_potencia_2, siguiente_potencia_2, dim(imagen)[3])) 
  
  imagen_cuadrada[1:n_filas, 1:n_columnas, ] <- imagen
  
  return(imagen_cuadrada)
}
>>>>>>> Alejandra





<<<<<<< HEAD
Suguiendo el mismo procedimiento, redimensionamos las imágenes y y
realizamos las transformadas wavelet y el thresholding con la función
procesar_imagen_wavelet.




images_recortadas <- lapply(images_for_imwd_cut_3, resize_imwd)



Warning: Assuming third dimension corresponds to colourWarning: Assuming third dimension corresponds to colourWarning: Assuming third dimension corresponds to colour



images_sin_ruidos<- lapply(images_recortadas, procesar_imagen_wavelet)
=======
Cargar imagenes en la función cuadrada
Aplicamos la función en 3 fotos




fotos_cuadradas <- lapply(imagen_noise, hacer_cuadrada_potencia_2)
>>>>>>> Alejandra




<<<<<<< HEAD
Visualizamos los resultados. Parece que el algoritmo funciona
correctamente para los tres tipos de ruido.




par(mfrow = c(3, 2), cex = 0.5)

display(Image(images_recortadas[[1]], colormode = 'Color'), method='r')
title('Imagen 2 con salt and pepper')




display(Image(images_sin_ruidos[[1]], colormode = 'Color'), method='r')
title('Imagen 2 sin ruido')





display(Image(images_recortadas[[2]], colormode = 'Color'), method='r')
title('Imagen 3 con ruido gamma')




display(Image(images_sin_ruidos[[2]], colormode = 'Color'), method='r')
title('Imagen 3 sin ruido')





display(Image(images_recortadas[[3]], colormode = 'Color'), method='r')
title('Imagen 3 con ruido uniforme')




display(Image(images_sin_ruidos[[3]], colormode = 'Color'), method='r')
title('Imagen 3 sin ruido')





=======
Imagen Original





>>>>>>> Alejandra




<<<<<<< HEAD
Redimensionamos las imágenes obtenidas a su tamaño original para
poder compararlas con estas. Vemos que en este caso la eliminación de
ruido ha sido muy buena y no hay apenas distorsión ni suavizado de los
bordes.




images_sin_ruido <- Map(resize_imwd_to_original, images_sin_ruidos, c("2", "3", "3") )








par(mfrow = c(3, 2), cex = 0.5)

display(Image(images[[2]], colormode = 'Color'), method='r')
title('Imagen 2 original')




display(Image(images_sin_ruido[[1]], colormode = 'Color'), method='r')
title('Imagen 2 sin ruido')





display(Image(images[[3]], colormode = 'Color'), method='r')
title('Imagen 3 original')




display(Image(images_sin_ruido[[2]], colormode = 'Color'), method='r')
title('Imagen 3 sin ruido')





display(Image(images[[3]], colormode = 'Color'), method='r')
title('Imagen 3 original')




display(Image(images_sin_ruido[[3]], colormode = 'Color'), method='r')
title('Imagen 3 sin ruido')











rm(images_recortadas, images_sin_ruido, images_sin_ruidos)




=======



[1] "1 Noise: gaussian"
[1] "2 Noise: gamma"







[1] "3 Noise: unif"













Cuando uso fotos de mayor calidad solo acepta hasta 2, con 3 aparece
esto: Error: vector memory limit of 16.0 Gb reached, see
mem.maxVSize()
El umbral “universal”, propuesta por Donoho y Johnstone, es un método
utilizado en la eliminación de ruido en señales e imágenes mediante
transformadas de wavelet. Esta estrategia calcula el umbral aplicado a
los coeficientes de wavelet en función del tamaño de la señal y una
estimación del nivel de ruido. La fórmula del umbral “universal” es
\[ \sigma \sqrt{2 \log n}\] donde
\(\sigma\) es una estimación del ruido
y n es el número de muestras o elementos de la señal.




Imagenes sin ruido



Error in slot(prueba, ".Data")[1:1600, 1:1066, , drop = FALSE] : 
  subscript out of bounds



Ruido gaussiano




Imagen con ruido Gaussiano original

Imagen Cortada sin filtro
Ajusta el contraste de la imagen, incrementando la diferencia
entre los píxeles claros y oscuros.
Imagen con filtro paso alto














4.3 Función denoise…
>>>>>>> Alejandra

8 Conclusiones

<<<<<<< HEAD ======= >>>>>>> Alejandra